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Abstract 

Context. The formation of vortices in accretion disks is of high interest in various astrophysical contexts, in particular 
for planet formation or in the disks of compact objects. But despite numerous attempts it has thus far not been possible 
to produce strong vortices in fully three-dimensional simulations of disks. 

Aims. The aim of this paper is to present the first 3D simulation of a strong vortex, established across the vertically 
stratified structure of a disk by the Rossby Wave Instability. 

Methods. Using the Versatile Advection Code (VAC), we set up a fully 3D cylindrical stratified disk potentially prone 
to the Rossby Wave Instability. 

Results. The simulation confirms the basic expectations obtained from previous 2D analytic and numerical works. The 
simulation exhibits a strong vortex that grows rapidly and saturates at a finite amplitude. On the other hand the third 
dimension shows unexpected additional behaviours that could be of strong importance in the astrophysical roles that 
such vortices can play. 

Key words. Accretion, accretion disks - Planetary systems: protoplanetary disks - Hydrodynamics - Instabilities - 
Methods: numerical 



1. Introduction 

Vortices have been actively searched for in accretion disk 
theory and numerical simulations because of their multiple 
astrophysical interests. In particular, in the case of proto- 
planetary disks, the presen ce of vortices has been invoked 
(Barge & Sommeria |1995[ ) to alleviate the problem of the 
too long time scale needed for the growth of grains to plan- 
etesimals. The streaming of the gas inside a vortex could 
speed up this growth by concentrating dust grains in its 
centre, by what is sometimes called the 'tea leaf effect 
whereby tea leaves accumulate in the center of a cup when 
it is stirred. This mechanism has been studied numerically 
( |Johansen et al.||2004| |Lyra et al.| |2008| |2009|) but only in 



the 2D case of an infinitely thin disk. We note that an al 
ternative mechanism for trapping the dust grains based on 
high p ressure region has been proposed by |Johansen et ah] 
p007 ). 



On the other hand, despite multiple attempts to gener- 
ate realistic vortices in fully cylindrical 3D numerical simu- 
lations of disk s, even the most advanced ones (Barranco & 
Marcus 2005 1 have only found o ff-midplane vortices. This 
may be explained ([Tagger 2001) by the fact that due to 
both differential rotation and differential vorticity (which 
is absent in the Shearing Box model used in many simu- 
lations) vortices are very rapidly sheared away and could 
survive only by being continuously re-generated. 

In this paper we present the first 3D numerical simula- 
tions of the Rossby Wave Instability (RWI), which is pre- 



cisely a mechanism to do this: in certain conditions this in- 
stability generates Rossby vortices that grow exponentially 
with time. The existence of this instability in protop lane- 



tary disks was first mention ed in Lovelace et al 
it has been suggested later (Varniere & Tagger 



( 1999) and 
2006) that 



at the edges of the 'Dead Zone' of protoplanetary disks, 
where the gas is not ionized and thus should not be tur- 
bulent, the conditions are such that the RWI can grow; it 
should thus both help accretion to proceed across the Dead 
Zone, and help planets to grow there. The simulations we 
present confirm the basic properties of the instability, but 
they show new features in its 3D structure, whose effects 
on accretion and the growth of planetesimals will need to 
be understood. 

The paper is organized as follows: after recalling the 
basic properties of the RWI we will present the physical 
setup, then the numerical one. Section 4 will present the 
results, and section 5 a discussion and conclusion. 



2. The Rossby Wave Instability 

The RWI has been found and discussed in different contexts 
of differentially rot ating disks. A discussion o f this history 
has been given by Varniere & Tagger (2006), an d we will 
not repeat it here. It can exist in galact ic disks (Lovelace 
|fc Hohlfeld||1978| |Sellwood fc Kahn|l99l|) , as well as in ac- 
cretion disks jPapaloizou &: Pringle| Q1985 

200ip ) 



Lovelace et al. 

|1999J; see also lLi et al.| ( |2000||2001| )TTtcan be seen as the 
form that the Kelvin-Helmoltz instability takes in differen- 
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Figure 1. Schematic view of the Rossby Wave Instability. 
Rossby waves are generated in the region of the density 
extremum, and density waves are emitted away from this 
region due to differential vorticity, which couples these two 
waves. The Rossby waves have their corotation radius at 
the density maximum 



tially rotating disks, and it has a similar instability crite- 
rion: the existence of an extremum of a quantity C related 
to the vorticity of the equilibrium flow. In an unmagnetiz ed 
thin disk this quantity can be written as ( Li et al. 2000 1 : 



C 



Eft p 

2^? £7 



P 



2(V x v) z Et 



(1) 



where E is the surface mass density, v is the velocity of 
the fluid, 7 the adiabatic index, f2 the rotation frequency 
and k 2 — + 2f2fi' the epicyclic frequency squared (so 
that K 2 /2f2 is the vorticity). Here the prime notes a radial 
derivative. 

Two possibilities that can result in an extremum in 
C have been investigated: the first one occurs near the 
marginally stable orbit around a compact object, where 
relativistic effects create a maximum of k 2 /2f2. The growth 
rate of the RWI is strongly in creased by a poloidal mag - 



netic field threading the disk (Tagger fe Varniere 2006) 
but decreased by a toroidal one ( Yu & Li 2009p . This has 
been proposed as an explanation for the high fr equency 
quasi-periodic oscillation (QPO) of microquasars flTag ger| 



fe Varniere 2006; L ai fe Tsang||2009| |Tsang fe Lai| [2009 

The second possibility is to have an ex tremum of the 
surface density. The model of Varniere & Tagger ( 2006 ) for 
protoplanetary disks relies on the fact that extrema of E 
should occur at the edges of the Dead Zo ne of these disks, 
so tha t the RWI should be unsta ble there. |Tagger fc Melia 



(2006) and Falanga et al. (2007) have also used the RWI 



in an explanation for the quasi-periodicity that may have 
been observed during the flares of SgrA*. 



The RWI has been stu died both analyti cally (Lovelace 
et al.||1999||Li et al ||2000[ ) and numerically ( |Li et al poTTf] 



V 



armcrc 



& Ta- 



gger 



2006 ) . It is formed by Rossby waves 



trapped in the extremum of C, as shown in figure [T] In the 
/3-plane approximation differential rotation is neglected and 
the dispersion relation of Rossby waves is given by 



ma 



k 2 



(2) 



where k is the adimensional radial wavenumber, m the az- 
imuthal wavenumber, and a the vorticity gradient. In an 



accretion disk, differential vorticity and differential rota- 
tion couple them to spiral waves, emitted on both sides of 
the extremum region. The wave dis persion relation calcu- 
late d from the equ ations obtained by Lovelace et al. ( 1999 ) 
andlTaggerl (|2001i) is 



c 2 (fc 2 + m?) + r 2 K 2 



(3) 



where w = uj — mtt, c s is the sound speed and r the radius. 

Rossby waves have their corotation radius (where their 
phase velocity equals the rotation velocity of the gas) at 
that extremum. The standing wave pattern they form ap- 
pears as a vortex located in the region of the extremum, 
and grows exponentially. The waves have a positive flux of 
energy and angular momentum beyond that radius, and a 
negative flux inside it, so that (as will be checked in the 
simulation) the pattern can grow as it causes the gas inside 
corotation to lose angular momentum and accrete, while 
the gas beyond corotation gains that angular momentum 
and moves outward. 

As a result the instability tends to destroy the ex tremum 
of C, as seen for instance in Tagger fc Melia ( 2006). In the 
contex t of the mechanism proposed by |Varniere fc Tagger| 
(2006) for protoplanetary disks, we expect the instability 



to saturate at the amplitude where this is balanced by the 
continuous regeneration of the extremum of density by the 
gas accreting from the outer disk and piling up at the edge 
of the Dead Zone (and oppositely at the inner edge of the 
Dead Zone). 

However these studies have all been done in the approx- 
imation of an infinitely thin disk, because of the complexity 
of a full 3D analytical study, or of the numerical resources 
needed for 3D simulations. Such a study is however highly 
desirable, both to consider the full complexity of the gas 
(and grains) flow, and to validate the thin disk approxi- 
mation: this approximation was introduced (and its limits 
defined) by Goldreich & Lynden-Bell ( 1965a|b[ ), but this 
does not apply as such to Rossby perturbations. We will 
see farther in this paper that indeed 3D effects bring in sig- 
nificant and potentially important qualitative differences. 



3. Gaseous accretion disk 

3.1. Hydrodynamics equations 

We work in cylindrical coordinates (r, </>, z) with the 3D 
Euler equation: 



d tP + V.{vp) = 
d t (pv) + V.(vpv) + Vp = -pV&c 



(4) 
(5) 



where p is the mass density of the fluid, and v its velocity, 
and p the pressure. $g = — GAf*/(r 2 + z 2 ) 1 / 2 is the grav- 
ity potential of the central object with G the gravitational 
constant and the mass of the central object. We con- 
sider a barotropic flow, i. e. the entropy S is constant in the 
entire system: 



p = Sp 1 



(0) 



with the adiabatic index 7 = 5/3. The sound speed is given 
by c 2 = 7p/p = S^p 7-1 and the temperature by T = p/p = 

Sp-1- 1 . 
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Figure 2. Isocontours of the initial density, with a bump 
at tb = 3r^. The density is normalized to the density in the 
midplane at the inner edge. 
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Figure 3. The initial azimuthal velocity of the midplane 
gas (continuous line) is compared to keplerian rotation 
(dots). A zoom on the bump region is shown in the up- 
per right corner. 

3.2. Initial conditions 

We choose as initial equilibrium a density profile decreasing 
radially as r" 1 / 2 , to which is added a density bump: 



p ml (r,z = 0) 



= 1 



X exp ■ 



(r/r.i - r B /rif 



2a 2 



1 



(7) 



The density is normalized by p m , the midplane density 
at the inner boundary of the simulation r^, and tb is the po- 
sition of the bump. For the parameters \ an d c, which con- 
trol the amplitude and the width of the density bump, we 
take the values x = 0-4 and a = 0.1 respectively that gives 
a bum p similar to the one obtained by |Varniere fc Tagger| 
(2006) in their simulation of the dead zone. These values 
allow the instability to grow relatively fast, decreasing the 



Figure 4. Plot of the logarithm of the square of the 
epicyclic frequency k 2 (averaged over z) as a function of 
the radius. 



numerical load, while the effect of the pressure gradient on 
the equilibrium rotation velocity leaves us a margin with 
respect to the Rayleigh criterion n 2 > (see figure Q. 

The bump is centered at r B = 3r^, far enough from 
the inner edge of the disk to avoid a strong effect of the 
boundary condition there. The vertical density profile is 
initially at hydrostatic equilibrium, giving an aspect ratio 
of the order of 10 _1 . 



p ini (r, z )=p(r,z = 0) 



1 - 



7' 



i 



~/Sp(r,z = O)''- 1 



(8) 



\/r 2 + z 2 



with and S = 10" 3 r 2 ^ 



Finally, we use initially 
a 'floor' density p m i„ = 10~ 2 p m in order to avoid getting 
too low densities in the corona. Figure [2] shows the resulting 
isodensity contours in a vertical cut of the disk. 

The initial velocity field is purely toroidal, u* m = vl m = 



and 



has been chosen for the disk to be in radial 



equilibrium in a Newtonian potential: 



" 2 GM. M 



[r 2 + z 2 ) 3 / 2 pn 



(9) 



Figure [3] shows the deviation of the disk azimuthal velocity 
from a pure keplerian rotation. The pressure term in the 
last equation induces only a weak variation. 

4. Numerical setup 

4.1. Numerical code and scheme 

The numerical simulations presented h ere use the V ersatile 
Advection Code (VAC) developped by |Toth| ( [T996| . In the 
version we use the code solves the 3D hydrodynamics equa- 
tions for an isentropic flow. We use VAC with the total vari- 
ation diminishing monotonic upstream-centered scheme for 
conservation laws (TVD-MUSCL), a Roe Riemann solver, 
a Hancock predictor step and a Woodward limiter ( Colella 
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& Woodward 198 4). T he TVD-MUSCL scheme detailed in 
T6th& Odstrcil (1996) is one of the less dissipative schemes 
included in the VAC code, which is useful in order to ob- 
serve the full development of the instability and properly 
characterize its saturation. 



4.2. Grid and numerical resolution 

The grid is cylindrical (r,4>,z) with a resolution of 154 x 68 x 
68, including 2 ghosts cells at each boundary in order to 
impose the boundary conditions. The grid is non uniform 
with a streching factor of 20 in the radial direction, and 
10 in the vertical one. This achieves a higher resolution 
in the regions of physical or numerical interest. The mesh 
extends from (ri,0,0) to (6r.;, r^, 2tt) where r and z have 
been normalized by the disk inner edge radius rj . The radial 
and vertical extensions are large enough for the different 
waves of the RWI to develop without unwanted effects from 
the boundary conditions. 

The minimum value of the mass density is fixed to 
10~ 6 p m . With the initial conditions presented before this 
configuration gives a bump with approximatively 60 cells 
radially and 40 vertically. 

4.3. Boundary conditions 

The boundary conditions are imposed using the two ghost 
cells added at each boundary of the mesh, allowing to com- 
pute the derivatives in all the cells. Since the simulation is 
restricted to a small part of the disk, we have chosen con- 
tinuous boundary conditions in the radial direction, mean- 
ing that the values of the physical quantities are copied 
from the first cell of the computational domain in the ghost 
zones. This boundary condition is thus partially transpar- 
ent. In the azimuthal direction the boundary condition is 
periodic. Only the upper part of the disk is simulated since 
the vortices we are interested in are even in z. In section [B~2l 
we discuss a test run without this condition, which confirms 
the validity of this choice. 

4.4. Disk equilibrium and stability 

As mentioned in section 13.21 we start the simulation with 
a low but finite density in the corona. On the other hand, 
we need to start from a true disk equilibrium in order to 
avoid a rapid relaxation that would dominate the evolution 
of the system. The main difficulty is the constant density 
we choose to start from in the corona, and the transition to 
the power law vertical profile in the disk; this introduces a 
jump in the vertical derivative of the density. We thus ap- 
ply to this density profile a smoothing over n smoo t/j vertical 
grid points (with n smoo th — 4 here) in order to enforce a 
smooth transition. All this means that we are not at exact 
hydrostatic equilibrium in the gravity field. A different but 
related difficulty arises from the discretization of the grid 
and the use of slope limiters, which introduces residual spu- 
rious forces for each cell but a zero mean value on the whole 
grid. Although they are small (at the level of roundoff error) 
they are systematic (i.e. always act in the same direction) 
and become important in low-density regions, especially at 
the disk-corona transition. 

We thus choose to slightly modify the gravitational po- 
tential, so that the initial pressure profile given in section 
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Figure 5. Plot of the normalized density in the midplane of 
the disk at t — 176/ flx(ri)- One can identify the position 
of the crescent-shaped vortex at the bump radius (7' = 3r^) 
and the spiral arms on both side of the bump region. 

|3.2| is really at equilibrium. We apply the following proce- 
dure: a first time step is done to calculate the difference 
between the physical and numerical equilibrium. This al- 
lows us to compute the source term pV$c in the Euler 
equations that enforces numerical equilibrium: this source 
term should be equal to the gravity term if we were in exact 
equilibrium, and the difference arises from these spurious 
forces. We enforce equilibrium by subsequently consider- 
ing this constant modified source term. We will a posteri- 
ori check that this does not substantially affect our result. 
Moreover we have checked that with the whole simulation 
box inside the disk (without corona) this modified source 
term is not needed, but wave reflection at the upper bound- 
ary makes it impossible to study properly the physics of the 
instability. 

After several time steps we add to these initial condi- 
tions low level perturbations so as to provide seeds for the 
unstable modes that can develop if the disk is unstable. We 
find more convenient to do this with perturbations of low 
azimuthal wavenumber, since (as confirmed in the simula- 
tion) they are the only ones expected to be unstable. We 
thus apply an initial radial velocity perturbation: 

v rp = v r + esin(27r(r - 1.2)/.8) [sin(<?3) + sm(2<£) (10) 

+ sin(3</>) + sin(50)] 

with the amplitude of the perturbation e of the order of 
10- 7 . 

5. Results 

In this section we describe the results obtained with this 
simulation. We first describe the general properties of the 
instability, which allow us to identify it and to study its 
structure and evolution. We then focus on the velocity 
stream to understand the flow pattern induced by the in- 
stability in the disk. 
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Figure 6. The compressional and rotational parts of the total flow (left) and its non-axisymmetric part (right) at 
t = 176/fiK-(ri), showing the two waves that form the instability. Top: the amplitude of (V x v) z traces the vorticity, 
i.e. the Rossby vortex centered at the extremum of density. Bottom: the amplitude of V.t; traces the compressional 
waves (spiral density waves) that develop away from that region. In these figures, the whole disk is represented with the 
azimuthal coordinate on the vertical axis and the radius on the horizontal one. Note that, although the top-right panel 
shows both a cyclonic and an anticyclonic vortex (as expected for an m — 1 mode), once combined with the bulk flow 
only the anticyclonic vortex appears in the top-left panel. 



5.1. General study of the instability 

Since the instability has only been studied in 2D, we will 
first base the discussion on what we obtain in the midplane 
of the 3D simulation. As expected the simulation shows af- 
ter a few rotation times a coherent perturbed pattern form- 
ing in the region of the density bump. Figure [5] shows the 
normalized density in the midplane. The non-axisymmetric 
pattern corresponds to the RWI, extending spiral arms on 
both sides of the bump. 

As explained in section[2] the structure of the instability 
is formed by Rossby waves on both sides of the extremum 
of C (here, the density bump at r = 3rj), resulting in a 
standing vortex there, and spiral density waves emitted on 
both sides of that region. We can visualize this by plotting 
the curl and divergence of the flow, showing respectively 
its rotational (Rossby waves) and compressional (density 
waves) parts. Figure [6] shows that the expected patterns do 
indeed appear. 

The different stages of the development of the instabil- 
ity can be seen on figure [7J which presents on a logarithmic 
scale the amplitude of the density perturbation as a func- 
tion of time. After a short stage where the non pertinent 
perturbations decrease, the unstable mode enters the linear 
stage, i.e. the perturbation grows exponentially. We find 



a growth rate of .39^x(r = ri)/Q.K{r = 3^). This is com- 
parable wit h the growth rates obtai ned from 2D theory ( Li 
et al.|2000 ) or numerical simulation ( Tagger & Melia 2006 ) . 
The different physical and numerical setups (including the 
equation of state and the boundary conditions) prevent a 
more detailed comparison. The local dispersion equation 
(eq. [3]) does not allow to estimate the growth rate of the 
global modes. This would be possible anyway only with a 
fully 3D linear computation, which would require a very 
important effort. The perfect exponential growth we ob- 
tain in the present non-linear simulation over four decades 
in amplitude, its coherence with the 2D results, and its in- 
dependence on initial conditions show that we do capture 
properly here the linear phase of the instability. This stage 
lasts about 10 keplerian orbits at the bump radius, and fi- 
nally the instability reaches saturation. The entire progres- 
sion of the instability is thus captured in this simulation. 

Figure [9] shows that the saturation is not due to the 
flattening of the initial density bump and thus of the ex- 
tremum of £, as in the 2D simulations of Tagger & Melia 
(2006), so that non-linearities must be suspected. In order 
to check this we have run a new simulation, taking as ini- 
tial radial density profile the average one at the end of the 
present simulation (but again with only weak velocity per- 
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Figure 7. Instability growth: the maximum value of the 
non axisymmetric part of the density is plotted as a func- 
tion of time, with a logarithmic scale for the vertical axis. 
We have here added a second time axis (lower) that shows 
the number of keplerian orbit at the bump radius in or- 
der to estimate the growth rate in the same units as 
Li et al. ( |2000[ ). The equation of the continuous line is 
— 4.1 + 0. 39n7^(ri)/r2 x (3ri). The graph thus shows the lin- 
ear stage and the saturation after about 10 keplerian times 
at r = 2>ri. 



turbation, as at its startup), and observed the RWI growing 
again. 

The azimuthal density profile at the bump radius is 
shown in figure 10. It reveals a departure from a pure sinu- 
soid, corresponding to the presence of m > 1 modes, whose 
growth is shown in figure 7. Since they do not show the ex- 
ponential growth of the m — 1 but they show modes that 
are not present in the initial perturbation, we analyze them 
as non-linearly generated harmonics of that mode. Figure 8 
shows that their effect remains relatively weak, so that the 
saturation of the mode cannot be ascribed to them. 

On the other hand, as shown in figure [l3j small-scale 
structures appear in the vertical flow, and become stronger 
when the mode approaches saturation. They cascade to the 
smallest scale of the simulation, making us believe that dis- 
sipation in these structures is the main cause of the satura- 
tion we observe. Assessing their role will require both im- 
proved numerics, to deal with these structures, and physical 
understanding since one can expect that vertically propa- 
gating sound waves generated within the disk will be sub- 
ject to wavebreaking when they reach the low-density, low- 
sound speed corona. This will be considered in future work. 

The flattening of the bump (figure [9]) is explained by 
the accretion rate, shown on figure |11| it is positive inward 
from the density extremum and negative beyond it. This 
is expected from a mode with its corotation at that radius 
which grows by exchanging angular momentum between the 
inner and outer region. As in 2D simulations, the transition 
between positive and negative accretion is so sharp that it 
is actually the best diagnostic of the corotation radius and 
thus of the mode frequency. 

Finally, considering the possibility of other unstable 
modes, we show in figure [12] the time evolution of the az- 
imuthal Fourier components, m = to 3. Although per- 
turbations do appear at m = 2 and 3, the fact that they 
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Figure 9. Comparison between the density integrated over 
z and (f> (left) and of the RWI criterion C (right) at t = in 
solid line and t — 176/Qjc( r i) m dashed line. One can see 
that the extremum of C at r = 3r^ has decreased but not 
disappeared when the instability saturates. 
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Figure 10. The azimuthal surface density profile at r = 3r^ 
at different times with the same normalization. The line is 
at t = 32/f2if(rj), the departure from the sinusoid comes 
from the initial perturbations that is a sum of sinus. The 
dots is at t = 104/Oif(ri), the profil is closer form a si- 
nusoid, the instability is in the linear phase, the nonperti- 
nent modes (m>l) are negligeable. The dashed line is at 
t = 176/fi#-(rj), the departure from a sinusoid is believed 
to be a direct consequence of non-linearities. 



evolve together with the m = 1 indicates that they are 
probably harmonics of that mode and not independently 
growing ones. The most linearly unstable m depends on 
the local conditions. We have run 2D simulations with the 
same density profile and find in that case a dominant m = 2. 
We will see below that the 3D simulation does show new 
features in the structure of the mode, which seem to con- 
tribute to its linear growth. They might be at the origin of 
this difference. 
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Figure 8. The evolution of the amplitude of the ten first azimuthal modes in the same coordinates as figure |7| The 
amplitude of each mode is calculated with an azimuthal Fourier transform, pi being the amplitude of the i th mode. The 
development of modes that were not present in the initial perturbation proves the existence of non-linearities. 



5.2. Study of the velocity stream 

In this section we consider the flow pattern generated by the 
instability. Since the stream in the (r, (f>) plane is dominated 
by keplerian rotation, we consider the non axisymmetric 
part of the velocity, i.e. we study the disk in a sheared 
frame rotating as the mean flow at each radius. Figure 13 
does show in this plane the expected vortex structure, with 
both a cyclonic and anticyclonic feature since the m = 1 
structure generates both. As seen in figure [5] however, when 
combined with the keplerian flow only the anticyclonic com- 
ponent appears, forming a crescent-shaped feature, while 
the cyclonic one is at the x-point where the tips of the cres- 
cent join. 

On the other hand figure [13] also shows unexpected new 
features: first, the plot in the (0, z) plane shows an down- 
draft at the center of the anticyclonic vortex, and a updraft 
at the center of the cyclonic one. Also, plots in the (r, z) 
plane show that strong convection rolls form on both sides 
of the vortex, in the anticlockwise direction for the outer 
one and clockwise direction for the inner one. Streamlines, 
at the outer edges of these horizontal and vertical vor- 
tices, connect them smoothly in a manner reminiscent of 



orbits near separatrices in hamiltonian dynamics. The main 
streamlines are schematized on figure [T4] 

We have checked the robustness of these convection rolls 
by performing numerical simulations with the entire verti- 
cal structure of the disk, without assuming a vertical sym- 
metry, all the other parameters kept unchanged. The verti- 
cal symmetry was broken by the initial perturbation added 
on the radial velocity, but resulted in similar flow patterns. 
Furthermore, when we applied only antisymmetric initial 
perturbations, we found that this delayed the development 
of the instability without changing its vertical structure. 
We conclude that the delay is due to the difficulty for anti- 
symmetric initial perturbations to seed the unstable mode, 
and that the convection rolls are part of the linear struc- 
ture of the instability. To confirm this diagnostic, we also 
checked that the maximal vertical velocity involved in these 
convection rolls, follows the linear growth of the instability 
as in figure [7] 

This strong vertical structure of the instability and of 
the vortex pattern could not be expected from the 2D re- 
sults. We believe that the vertical vortices are resonantly ex- 
cited where the local Doppler-shifted frequency, uj — mfl(r), 
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Figure 12. Temporal evolutions of the first four azimuthal modes. A different color bar is used for m=0 which is 
approximately ten times stronger than the others. 



accretion 




Figure 11. The accretion rate defined as 

— J rpv r dzd<t>/ J r p m c s dzd4> as a function of radius, Figure 14. Schematic view of the circulation, showing the 
at the end of the linear phase (t = 176/CIk (fi)) horizontal and vertical vortices in a sheared frame rotating 

at the local mean angular velocity. 



is equal to the mean vertical oscillation frequency. The rela- 
tion between this and the Brunt- Vaisala frequency, depend- 
ing on the equation of state of the gas, will be discussed else- 



where. This would thus be the equivalent, in a fluid disk, of 
the 'peanuts' structure found in observations and numeri- 
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cal simulations of barred spiral galaxies ([Combes fc Sanders 
1981| ICombes et al.|[l990{ |Athanassoula||2008l ). This is beT 
lieved (although an explanation based on the firehose insta- 
bility has also been proposed), to result from the resonant 
excitation of the vertical motion of stars where u) — mf2(r) 
equals their vertical frequency of motion. 



6. Discussion and outlooks 

We have presented here the first fully three-dimensional, 
cylindrical numerical simulation showing a strong and per- 
sistent vortex in a differentially rotating disk. Its growth 
results from the Rossby Wave Instability, which has been 
invoked in various astrophysical contexts such as in proto- 
planetary disks or in the accretion disk of compact objects. 

The simulation essentially confirms the main properties 
of the instability, which forms a vortex extending across 
the vertically stratified structure of the disk. The vortex is 
localized at the extremum of vortensity and very efficiently 
acts to destroy this extremum, by permitting an exchange 
of angular momentum between regions located on either 
side of it. 

In ou r simulation we find no sign of the ellipt ical in- 



stability flKerswell||2002} |Lesur k. Papaloizou||2009| ) which 
has been claimed to destroy vortices. We note however that 
these works are done in the shearing sheet approximation 
which neglects the equilibrium vorticity gradient and thus 



can consider local small-scale modes whereas the RWI is a 
large-scale global instability. The s ame remark applies to 
the work of [Latter & Balbusl <|2009b. 



The most important new result found in our 3D simu- 
lation is the occurrence of strong vertical convection rolls, 
excited on both sides of the main, horizontal, vortex. Figure 
13] shows that the vertical velocity in these rolls is compa- 
rable to the radial velocity involved in the Rossby vortex, 
while their growth shows that they are an inherent part of 
the structure of the mode, tightly connecting the horizontal 
and vertical circulations. 

Future work should allow us to assess the importance 
of this convection in the astrophysical contexts where the 
RWI may act. In particular we expect them to play a very 
important role in the concentration of dust grains in proto- 
planetary disks, in addition to the mechanism of Barge fc] 



Sommeria ( 1995 ), which will accumulate them at the center 
of the horizontal vortex. Indeed figure [13] also shows that, in 
the midplane of the disk, the flow lines rapidly spiral inward 
in the cyclonic vortex, and outward in the anticyclonic one. 
This goes together with the updraft at the cyclonic vortex 
and downdraft at the anticyclonic one. We note that dust 
grains will thus be rapidly transported with the gas, at low 
z, toward the center of the cyclonic vortex but that their 
weight should make it very difficult for them to be dragged 
in the updraft. The 3D structure of the vortex presented 
here will then have a direct impact on the accumulation of 
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grains that will happen on an even lower timescale than in a 
2D vortex. We also note howeve r that, counter to both intu - 
ition and the effect discussed by Barge & Sommeria ( 1995 1, 
this would tend to accumulate the grains at the cyclonic 
rather than anticyclonic vortex. 

In order to explore this mechanism, which could be of 
primordial importance for the growth of planctesimals and 
planet formation, work should proceed in two directions: 
performing simulations using varied density profiles, includ- 
ing the one expect ed at the edge of the Dea d Zone of pro- 
toplanetary disks ( Varniere fc Taggerj|2006[ ), and following 
the motion and growth of dust grains in the resulting verti- 
cal flows. An analytical description of the vertical structure 
observed in the simulation is also needed. 

Another direction for future work concerns the initial 
goal of this work, which was to analyze the 3D struc- 
ture of MHD instabilities in magnetized disks. In accre- 
tion disks threaded by a vertical magnetic field, it has been 
sho wn that an Accretion- E jection Instability (AEI) can oc- 



cur (Tagger & Pellat 



1999), and explain the lo w-frequency 



Quasi-Periodic Oscillation of X-r ay binaries ( | Rodriguez 
et al. 2002 Varniere et al. 2002 1. The instability grows 



by coupling magnetically-driven spiral density waves and a 



Rossb y vortex. It has also been shown (Varniere & Tagger 
2002) that this vortex can re-emit vertically, as Alfven 



waves propagating to the corona along magnetic field lines, 
a substantial fraction of the accretion energy and angu- 
lar momentum; it was also suggested that this could be 
a source for accretion-driven winds and jets. In this con- 
text we can expect that the convection rolls observed here 
would be replaced by resonantly excited slow magnetosonic 
waves. Since the main action of these waves is to move gas 
along magnetic field lines we can expect them to provide 
mass-load to the resulting jet, much as in MHD steady-state 



models of jets ( Ferreira & Pelletier 



1995 Casse & Ferreira 



2000) mass-loading occurs at the slow magnetosonic point 
and acceleration occurs in the region of the Alfvenic point. 
In order to study this prospect a fully 3D simulation that 
can handle both the geometry (a disk threaded by open 
poloidal magnetic field lines) and the magnitude (of the or- 
der of equipartition with the gas pressure) of the relevant 
magnetic field needs to be developed. 
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